#######################################
## Replication code for
## "Social Contact and Outgroup Bias: World Cup Matches and Changes in Outgroup Bias"
## Jong Hee Park and Anakiz Elif Senturk
#######################################

################################################
## helper packages and codes
################################################
library(tidyverse)
library(haven)
library(sjlabelled)
library(labelled) 
library(surveytoolbox) 
library(dplyr)
library(list)
library(dotwhisker)
library(gtsummary)
library(gridExtra)
library(list)
library("cobalt")
library("MatchIt")



model.maker <- function(input){
    data <- rbind(cbind(coef(input)[2], sqrt(diag(vcov(input))[2])),
                     cbind(coef(input)[1], sqrt(diag(vcov(input))[1])))
    data.frame("term" = c("Control", "Outgroup"),
               "estimate" = data[, 1],
               "std.error" = data[, 2],
               "statistic" = data[, 1]/data[, 2],
               "p.value" = pt(data[, 1]/data[, 2], df = c(length(input$y), length(input$y)), lower.tail=FALSE), 
               "model" = c("Control", "Outgroup"))%>% as_tibble()
}

ict.table.maker <- function(input){
    res <- cbind(input$par.treat, input$se.treat)
    colnames(res) <- c("coef", "se")
    return(res)
}


################################################
## call data
################################################
df <- readRDS("Data/replication_data.rds")

################################################
## descriptive stat
## output is a messy latex code that needs editing
################################################
tbl_strata_ex1 <-
  df %>%
  dplyr::select(sex, residence7, partisanship, age.cat,
                political_ideology, 
                education, monthly_income_recode, gb) %>%
  gtsummary::tbl_strata(
                 strata = gb,
                 .tbl_fun =
                     ~ .x %>%
                         tbl_summary() %>% ## , missing = "no") %>% by = source
                         add_n(),
                 .header = "**{strata}**, N = {n}"
             )
tbl_strata_ex1%>%
    bold_labels()%>%
    as_gt() %>%
    gt::as_latex()


####################################################
## List experiment analysis
####################################################
df.A <- df %>%
    mutate(treat = ifelse(gb==1, 0, 1)) %>%
    ## treatment 1: south american
    mutate(y1 = ifelse(!is.na(y_listexp_control), y_listexp_control,
                ifelse(!is.na(y_stereotype_samerican), y_stereotype_samerican,
                       NA))) %>%
    mutate(treatment1 = ifelse(!is.na(y_listexp_control), 0,
                        ifelse(!is.na(y_stereotype_samerican), 1,
                               NA)))%>%
    ## treatment 2: african
    mutate(y2 = ifelse(!is.na(y_listexp_control), y_listexp_control,
                ifelse(!is.na(y_stereotype_african), y_stereotype_african,
                       NA))) %>%
    mutate(treatment2 = ifelse(!is.na(y_listexp_control), 0,
                        ifelse(!is.na(y_stereotype_african), 1,
                               NA)))%>%
    
    ## treatment 3: white
    mutate(y3 = ifelse(!is.na(y_listexp_control), y_listexp_control,
                ifelse(!is.na(y_stereotype_white), y_stereotype_white,
                       NA))) %>%
    mutate(treatment3 = ifelse(!is.na(y_listexp_control), 0,
                        ifelse(!is.na(y_stereotype_white), 1,
                               NA)))%>%
    
    ## treatment 1: outgroup bias
    mutate(y.binary = ifelse(!is.na(y_listexp_control), y_listexp_control,
               ifelse(!is.na(y_stereotype_samerican), y_stereotype_samerican,
               ifelse(!is.na(y_stereotype_african), y_stereotype_african,
               ifelse(!is.na(y_stereotype_white), y_stereotype_white,
                      NA))))) %>%
    mutate(treatment.binary = ifelse(!is.na(y_listexp_control), 0,
                       ifelse(!is.na(y_stereotype_samerican), 1,
                       ifelse(!is.na(y_stereotype_african), 1,
                       ifelse(!is.na(y_stereotype_white), 1,
                              NA)))))
 
################################################
## Calculate list experiment effect as difference in means
################################################
ict.results <- list()
ict.results[[1]] <- ictreg(y1 ~ 1, data = data.frame(df.A), 
                           treat = "treatment1", J=5, method = "lm")
ict.results[[2]] <- ictreg(y2 ~ 1, data = data.frame(df.A), 
                           treat = "treatment2", J=5, method = "lm")
ict.results[[3]] <- ictreg(y3 ~ 1, data = data.frame(df.A), 
                           treat = "treatment3", J=5, method = "lm")
ict.results[[4]] <- ictreg(y.binary ~ 1, data = data.frame(df.A), 
                           treat = "treatment.binary", J=5, method = "lm")


## R0
ict.results.R0 <- list()
ict.results.R0[[1]] <- ictreg(y1 ~ 1, data = data.frame(df.A %>% filter(gb==1)), 
                           treat = "treatment1", J=5, method = "lm")
ict.results.R0[[2]] <- ictreg(y2 ~ 1, data = data.frame(df.A %>% filter(gb==1)), 
                           treat = "treatment2", J=5, method = "lm")
ict.results.R0[[3]] <- ictreg(y3 ~ 1, data = data.frame(df.A %>% filter(gb==1)), 
                           treat = "treatment3", J=5, method = "lm")
ict.results.R0[[4]] <- ictreg(y.binary ~ 1, data = data.frame(df.A %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")

## R1
ict.results.R1 <- list()
ict.results.R1[[1]] <- ictreg(y1 ~ 1, data = data.frame(df.A %>% filter(gb==2)), 
                           treat = "treatment1", J=5, method = "lm")
ict.results.R1[[2]] <- ictreg(y2 ~ 1, data = data.frame(df.A %>% filter(gb==2)), 
                           treat = "treatment2", J=5, method = "lm")
ict.results.R1[[3]] <- ictreg(y3 ~ 1, data = data.frame(df.A %>% filter(gb==2)), 
                           treat = "treatment3", J=5, method = "lm")
ict.results.R1[[4]] <- ictreg(y.binary ~ 1, data = data.frame(df.A %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")

## R2
ict.results.R2 <- list()
ict.results.R2[[1]] <- ictreg(y1 ~ 1, data = data.frame(df.A %>% filter(gb==3)), 
                           treat = "treatment1", J=5, method = "lm")
ict.results.R2[[2]] <- ictreg(y2 ~ 1, data = data.frame(df.A %>% filter(gb==3)), 
                           treat = "treatment2", J=5, method = "lm")
ict.results.R2[[3]] <- ictreg(y3 ~ 1, data = data.frame(df.A %>% filter(gb==3)), 
                           treat = "treatment3", J=5, method = "lm")
ict.results.R2[[4]] <- ictreg(y.binary ~ 1, data = data.frame(df.A %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")

## R3
ict.results.R3 <- list()
ict.results.R3[[1]] <- ictreg(y1 ~ 1, data = data.frame(df.A %>% filter(gb==4)), 
                           treat = "treatment1", J=5, method = "lm")
ict.results.R3[[2]] <- ictreg(y2 ~ 1, data = data.frame(df.A %>% filter(gb==4)), 
                           treat = "treatment2", J=5, method = "lm")
ict.results.R3[[3]] <- ictreg(y3 ~ 1, data = data.frame(df.A %>% filter(gb==4)), 
                           treat = "treatment3", J=5, method = "lm")
ict.results.R3[[4]] <- ictreg(y.binary ~ 1, data = data.frame(df.A %>% filter(gb==4)), 
                           treat = "treatment.binary", J=5, method = "lm")

################################################
## all summary
################################################
ict.out <- rbind(cbind(coef(ict.results[[1]])[2], sqrt(diag(vcov(ict.results[[1]]))[2])),
                 cbind(coef(ict.results[[1]])[1], sqrt(diag(vcov(ict.results[[1]]))[1])),
                 cbind(coef(ict.results[[2]])[1], sqrt(diag(vcov(ict.results[[2]]))[1])),
                 cbind(coef(ict.results[[3]])[1], sqrt(diag(vcov(ict.results[[3]]))[1])),
                 cbind(coef(ict.results[[4]])[1], sqrt(diag(vcov(ict.results[[4]]))[1])))

ict_model <- data.frame("term" = c("Control", "South American", "African", "White", "Outgroup (all)"),
                       "estimate" = ict.out[, 1],
                       "std.error" = ict.out[, 2],
                       "statistic" = ict.out[, 1]/ict.out[, 2],
                       "p.value" = pt(ict.out[, 1]/ict.out[, 2],
                                      df = c(length(ict.results[[1]]$y),
                                             length(ict.results[[1]]$y),
                                             length(ict.results[[2]]$y),
                                             length(ict.results[[3]]$y),
                                             length(ict.results[[4]]$y)),
                                      lower.tail=FALSE), 
                       "model" = c("Control", "South American", "African", "White", "Outgroup"),
                       ict.out)%>% as_tibble()

################################################
## round by round summary
################################################
ict.out.R0 <- rbind(cbind(coef(ict.results.R0[[1]])[2], sqrt(diag(vcov(ict.results.R0[[1]]))[2])),
                    cbind(coef(ict.results.R0[[1]])[1], sqrt(diag(vcov(ict.results.R0[[1]]))[1])),
                    cbind(coef(ict.results.R0[[2]])[1], sqrt(diag(vcov(ict.results.R0[[2]]))[1])),
                    cbind(coef(ict.results.R0[[3]])[1], sqrt(diag(vcov(ict.results.R0[[3]]))[1])),
                    cbind(coef(ict.results.R0[[4]])[1], sqrt(diag(vcov(ict.results.R0[[4]]))[1])))

ict_model.R0 <- data.frame("term" = c("Control", "South American", "African", "White", "Outgroup"),
                           "estimate" = ict.out.R0[, 1],
                           "std.error" = ict.out.R0[, 2],
                           "statistic" = ict.out.R0[, 1]/ict.out.R0[, 2],
                           "p.value" = pt(ict.out.R0[, 1]/ict.out.R0[, 2],
                                          df = c(length(ict.results.R0[[1]]$y),
                                                 length(ict.results.R0[[1]]$y),
                                                 length(ict.results.R0[[2]]$y),
                                                 length(ict.results.R0[[3]]$y),
                                                 length(ict.results.R0[[4]]$y)),
                                          lower.tail=FALSE), 
                           "model" = c("Control", "South American", "African", "White",  "Outgroup"),
                           ict.out.R0)%>% as_tibble()


ict.out.R1 <- rbind(cbind(coef(ict.results.R1[[1]])[2], sqrt(diag(vcov(ict.results.R1[[1]]))[2])),
                    cbind(coef(ict.results.R1[[1]])[1], sqrt(diag(vcov(ict.results.R1[[1]]))[1])),
                    cbind(coef(ict.results.R1[[2]])[1], sqrt(diag(vcov(ict.results.R1[[2]]))[1])),
                    cbind(coef(ict.results.R1[[3]])[1], sqrt(diag(vcov(ict.results.R1[[3]]))[1])),
                    cbind(coef(ict.results.R1[[4]])[1], sqrt(diag(vcov(ict.results.R1[[4]]))[1])))

ict_model.R1 <- data.frame("term" = c("Control", "South American", "African", "White", "Outgroup"),
                           "estimate" = ict.out.R1[, 1],
                           "std.error" = ict.out.R1[, 2],
                           "statistic" = ict.out.R1[, 1]/ict.out.R1[, 2],
                           "p.value" = pt(ict.out.R1[, 1]/ict.out.R1[, 2],
                                          df = c(length(ict.results.R1[[1]]$y),
                                                 length(ict.results.R1[[1]]$y),
                                                 length(ict.results.R1[[2]]$y),
                                                 length(ict.results.R1[[3]]$y),
                                                 length(ict.results.R1[[4]]$y)),
                                          lower.tail=FALSE), 
                           "model" = c("Control", "South American", "African", "White", "Outgroup"),
                           ict.out.R1)%>% as_tibble()


ict.out.R2 <- rbind(cbind(coef(ict.results.R2[[1]])[2], sqrt(diag(vcov(ict.results.R2[[1]]))[2])),
                    cbind(coef(ict.results.R2[[1]])[1], sqrt(diag(vcov(ict.results.R2[[1]]))[1])),
                    cbind(coef(ict.results.R2[[2]])[1], sqrt(diag(vcov(ict.results.R2[[2]]))[1])),
                    cbind(coef(ict.results.R2[[3]])[1], sqrt(diag(vcov(ict.results.R2[[3]]))[1])),
                    cbind(coef(ict.results.R2[[4]])[1], sqrt(diag(vcov(ict.results.R2[[4]]))[1])))

ict_model.R2 <- data.frame("term" = c("Control", "South American", "African", "White", "Outgroup"),
                           "estimate" = ict.out.R2[, 1],
                           "std.error" = ict.out.R2[, 2],
                           "statistic" = ict.out.R2[, 1]/ict.out.R2[, 2],
                           "p.value" = pt(ict.out.R2[, 1]/ict.out.R2[, 2],
                                          df = c(length(ict.results.R2[[1]]$y),
                                                 length(ict.results.R2[[1]]$y),
                                                 length(ict.results.R2[[2]]$y),
                                                 length(ict.results.R2[[3]]$y),
                                                 length(ict.results.R2[[4]]$y)),
                                          lower.tail=FALSE), 
                           "model" = c("Control", "South American", "African", "White", "Outgroup"),
                           ict.out.R2)%>% as_tibble()


ict.out.R3 <- rbind(cbind(coef(ict.results.R3[[1]])[2], sqrt(diag(vcov(ict.results.R3[[1]]))[2])),
                    cbind(coef(ict.results.R3[[1]])[1], sqrt(diag(vcov(ict.results.R3[[1]]))[1])),
                    cbind(coef(ict.results.R3[[2]])[1], sqrt(diag(vcov(ict.results.R3[[2]]))[1])),
                    cbind(coef(ict.results.R3[[3]])[1], sqrt(diag(vcov(ict.results.R3[[3]]))[1])),
                    cbind(coef(ict.results.R3[[4]])[1], sqrt(diag(vcov(ict.results.R3[[4]]))[1])))

ict_model.R3 <- data.frame("term" = c("Control", "South American", "African", "White", "Outgroup"),
                           "estimate" = ict.out.R3[, 1],
                           "std.error" = ict.out.R3[, 2],
                           "statistic" = ict.out.R3[, 1]/ict.out.R3[, 2],
                           "p.value" = pt(ict.out.R3[, 1]/ict.out.R3[, 2],
                                          df = c(length(ict.results.R3[[1]]$y),
                                                 length(ict.results.R3[[1]]$y),
                                                 length(ict.results.R3[[2]]$y),
                                                 length(ict.results.R3[[3]]$y),
                                                 length(ict.results.R3[[4]]$y)),
                                          lower.tail=FALSE), 
                           "model" = c("Control", "South American", "African", "White", "Outgroup"),
                           ict.out.R3)%>% as_tibble()

ict.all <- bind_rows(ict_model %>% mutate(model="All"),
                     ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 

ict.out <- rbind(cbind(coef(ict.results[[1]])[2], sqrt(diag(vcov(ict.results[[1]]))[2])),
                 cbind(coef(ict.results[[1]])[1], sqrt(diag(vcov(ict.results[[1]]))[1])),
                 cbind(coef(ict.results[[2]])[1], sqrt(diag(vcov(ict.results[[2]]))[1])),
                 cbind(coef(ict.results[[3]])[1], sqrt(diag(vcov(ict.results[[3]]))[1])),
                 cbind(coef(ict.results[[4]])[1], sqrt(diag(vcov(ict.results[[4]]))[1])))

ict_model <- data.frame("term" = c("Control", "South American", "African", "White", "Outgroup"),
                        "estimate" = ict.out[, 1],
                        "std.error" = ict.out[, 2],
                        "statistic" = ict.out[, 1]/ict.out[, 2],
                        "p.value" = pt(ict.out[, 1]/ict.out[, 2],
                                       df = c(length(ict.results[[1]]$y),
                                              length(ict.results[[1]]$y),
                                              length(ict.results[[2]]$y),
                                              length(ict.results[[3]]$y),
                                              length(ict.results[[4]]$y)),
                                       lower.tail=FALSE), 
                        "model" = c("Control", "South American", "African", "White",  "Outgroup"),
                        ict.out)%>% as_tibble()
########################
## four hypotheses plots
########################
df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))


df.five.null <- df.five.alt <- df.five.slope <- df.five.valley <- df.five


df.five.null$estimate <- df.five.null$estimate[1]
df.five.alt$estimate[2:4] <- df.five.alt$estimate[3]
df.five.slope$estimate[2] <- df.five.slope$estimate[1]
df.five.slope$estimate[4] <- df.five.slope$estimate[1] - df.five.slope$estimate[3]
df.five.valley$estimate[2] <- df.five.valley$estimate[1]
df.five.valley$estimate[4] <- df.five.valley$estimate[3]



g1 <-  df.five.null %>% ## filter(term !="Treatment") %>%
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") + geom_point(size = 5, color="brown", alpha=0.5) +
    geom_point(data = df.five.null %>% filter(term == "Pre-game"), size = 5, color="brown", alpha=0.9) +
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    ggtitle("Null Hypothesis \n(No Update)") + ylim(-0.4, 0.7) +
    labs(y = "Outgroup Bias") +
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

g2 <-  df.five.alt %>% ## filter(term !="Treatment") %>%
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") + geom_point(size = 5, color="brown", alpha=0.5) +
    geom_point(data = df.five.null %>% filter(term == "Pre-game"), size = 5, color="brown", alpha=0.9) +
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    ggtitle("Attention Hypothesis \n(Update by Attention)") + ylim(-0.4, 0.7) +
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")
    
g3 <-  df.five.slope %>% ## filter(term !="Treatment") %>%
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") + geom_point(size = 5, color="brown", alpha=0.5) +
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five.null %>% filter(term == "Pre-game"), size = 5, color="brown", alpha=0.9) +
    ggtitle("Outcome Hypothesis \n(Update by Outcome)") + ylim(-0.4, 0.7) +
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")
   
g4 <-  df.five.valley %>% ## filter(term !="Treatment") %>%
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") + geom_point(size = 5, color="brown", alpha=0.5) +
    geom_point(data = df.five.null %>% filter(term == "Pre-game"), size = 5, color="brown", alpha=0.9) +
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    ggtitle("Unexpectedness Hypothesis \n(Update by Unexpected Outcome)") +
    ylim(-0.4, 0.7) +
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

pdf("four_hypothesis_plot.pdf", family="sans",
    width     = 14.5, height    = 3.5)
NetworkChange:::multiplot(g1, g2, g3, g4, cols=4)
dev.off()

########################
## actual outcomes
########################
df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))

g.main <-  df.actual %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five.null %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    ## geom_pointrange(aes(x=term, ymin=estimate.all-1.66*std.error.all, ymax=estimate.all+1.66*std.error.all),
    ##                color="brown", size=0.3, alpha=0.3) + 
    ylim(-0.4, 0.7) +
    labs(title = "(a) Main Result", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

pdf("result_plot.pdf", family="sans",
    width     = 8, height    = 6)
g.main
dev.off()

################################################
## actual outcomes with specific target countries only 
################################################
df.actual.target <- ict.all %>% filter(term=="Outgroup"& model=="Pre-game" |
                                       term=="South American"& model=="Round 1 (draw)" |
                                       term=="African"& model=="Round 2 (loss)" |
                                       term=="White"& model=="Round 3 (win)" )
df.actual.target$term[df.actual.target$term=="White"] <- "European"
df.actual.target$term <- factor(df.actual.target$term, levels=c("Outgroup", "African", "South American", "European"))
g.target <-  df.actual.target %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five.null %>% filter(term == "Pre-game") %>% mutate(term = "Outgroup"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    ylim(-0.4, 0.85) +
    labs(title = "(b) Target Outgroup Only", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

pdf("result_target_only_plot.pdf", family="sans",
    width     = 8, height    = 6)
g.target
dev.off()

## statistical results
sink("average_table.txt")
xtable(df.five)
sink()



#######################
## gender
##  female 1089   male 1080 
#######################
## male
df.sub <- df.A %>% filter(sex=="Male")
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)


ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 


df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))

df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))
g.ma <-  df.five %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Pre-game Outgroup Bias") + 
    labs(title = "(c-1) Gender: Male", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")



## statistical results
sink("male_table.txt")
xtable(df.five)
sink()



## female
df.sub <- df.A %>% filter(sex=="Female")
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)
ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 
df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))

df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))

g.f <-  df.five %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") +
    labs(title = "(c-2) Gender: Female", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

pdf("results_plot_gender.pdf", family="sans",
    width     = 12.5, height    = 5)
NetworkChange:::multiplot(g.ma, g.f, cols=2)
dev.off()


## statistical results
sink("female_table.txt")
xtable(df.five)
sink()

#######################
## political ideology
## 12 liberal
## 45 conservative
## 3 moderate
#######################
## liberal
df.sub <- df.A %>% filter(political_ideology<3)
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)
ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 
df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))

df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))
g.l <-  df.five %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    labs(title = "(d-1) Ideology: Liberal", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")


## statistical results
sink("liberal_table.txt")
xtable(df.five)
sink()


## Moderate
df.sub <- df.A %>% filter(political_ideology==3)
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)
ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 
df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))
df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))

g.m <-  df.five %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    labs(title = "(d-2) Ideology: Moderate", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")


## statistical results
sink("moderate_table.txt")
xtable(df.five)
sink()



## Conservative
df.sub <- df.A %>% filter(political_ideology>3)
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)
ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 
df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))
df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))

g.c <-  df.five %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    labs(title = "(d-3) Ideology: Conservative", y = "Outgroup Bias") +
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")


pdf("results_plot_ideology.pdf", family="sans",
    width     = 13.5, height    = 4.5)
NetworkChange:::multiplot(g.l, g.m, g.c, cols=3)
dev.off()

## N = 104 + 609 = 713, 806, 468 + 70 = 538


## statistical results
sink("conservative_table.txt")
xtable(df.five)
sink()


#################################################
## merged graph
#################################################
pdf("results_for_pnas.pdf", family="sans",
    width     = 12, height    = 10)
grid.arrange(g.main, g.target, g.ma, g.f, g.l, g.m, g.c,
             ncol = 6, 
             layout_matrix = rbind(c(1,1,1, 2, 2, 2),
                                   c(3,3,3, 4, 4, 4),
                                   c(5, 5, 6, 6, 7, 7)))
dev.off()

#######################
## Age
## 360 313 397 433 666
#######################
## young
df.sub <- df.A %>% filter(age_recode < 4)
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)


ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 


df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))

df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))
g.y <-  df.five %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    labs(title = "Experiment Results: Under 50s",
         subtitle = "Consistent with Unexpectedness Hypothesis", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

## Old
df.sub <- df.A %>% filter(age_recode > 3)
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)
ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 
df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))

df.actual <- df.five ## %>% filter(term !="Treatment")
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))

g.o <-  df.five %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    labs(title = "Experiment Results: Over 50s",
         subtitle = "Consistent with Unexpectedness Hypothesis", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

pdf("results_plot_age.pdf", family="sans",
    width     = 12.5, height    = 6)
NetworkChange:::multiplot(g.y, g.o, cols=2)
dev.off()

##############################################
## ceiling and floor effect
##############################################
df.A$age.young <- ifelse(df.A$age.cat < 4, 1, 0)
df.A$conservative <- ifelse(df.A$political_ideology>3, 1, 0)
df.A$male <- ifelse(df.A$sex =="Male", 1, 0)
df.A$college <- ifelse(df.A$education >2, 1, 0)
df.A$treat <- df.A$treatment.binary == 1
df.A1 <- df.A %>% dplyr::select(y.binary, male, college, age.young, conservative, treatment.binary, gb) %>% drop_na()
df.A1 <- data.frame(df.A1)


## Replicate Table 7 of Blair and Imai
## Fit standard design ML model
## Replicates Table 7 Columns 1-2 in Blair and Imai (2010)
noboundary.results <- ictreg(y.binary ~ male + college + age.young + conservative, 
                             data = df.A1, 
                             treat = "treatment.binary", J=5, method = "ml")

summary(noboundary.results)

## Replicates Table 7 Columns 3-4 in Blair and Imai (2010)
ceiling.results <- ictreg(y.binary ~ male + college + age.young + conservative, 
                        data = df.A1, 
                        treat = "treatment.binary", J=5, method = "ml", ## ffit.start = "nls",
                        ceiling = TRUE, ceiling.fit = "bayesglm",
                        ceiling.formula = ~ male + college + age.young + conservative)

summary(ceiling.results, boundary.proportions=TRUE)

# Fit standard design ML model with floor effects alone
# Replicates Table 7 Columns 5-6 in Blair and Imai (2010)

floor.results <- ictreg(y.binary ~ male + college + age.young + conservative, 
                        data = df.A1, 
                        treat = "treatment.binary", J=5, method = "ml", ## fit.start = "nls",
                        floor = TRUE, floor.fit = "bayesglm",
                        floor.formula = ~ male + college + age.young + conservative)

summary(floor.results, boundary.proportions=TRUE)

# Fit standard design ML model with floor and ceiling effects
# Replicates Table 7 Columns 7-8 in Blair and Imai (2010)

both.results <- ictreg(y.binary ~ male + college + age.young + conservative, 
                       data = df.A1, 
                       treat = "treatment.binary", J=5, method = "ml", ## ffit.start = "nls",
                       floor = TRUE, ceiling = TRUE, 
		       floor.fit = "bayesglm", ceiling.fit = "bayesglm",
		       floor.formula = ~ male + college + age.young + conservative,
		       ceiling.formula = ~ male + college + age.young + conservative)

summary(both.results, boundary.proportions=TRUE)


ict_model.R0 <- ict.table.maker(noboundary.results)
ict_model.R1 <- ict.table.maker(ceiling.results)
ict_model.R2 <- ict.table.maker(floor.results)
ict_model.R3 <- ict.table.maker(both.results)


ict.all <- bind_cols(ict_model.R0,
                     ict_model.R1,
                     ict_model.R2,
                     ict_model.R3) 
sink("ceiling_floor_result.txt")
print(xtable(ict.all, digits=3))
sink()


## ML estimates of the estimated population probability of lying separately for ceiling and floor effects.
## Ceiling liar pop. prop. est. = 3e-04 (s.e. = 0)
## Floor liar pop. prop. est. = 0.04288 (s.e. = 0.00075)
## Ceiling liar pop. prop. est. = 0.00036 (s.e. = 0)
## Floor liar pop. prop. est. = 0.04236 (s.e. = 0.00075)

################################################
## balance check across treatment conditions
################################################
## 1:1 NN matching w/ replacement on a logistic regression PS
m.out.all <- matchit(treatment.binary ~ male + college + age.young + conservative, data = df.A1,
                     replace = TRUE)

list.bal <- list()
list.bal[[1]] <- summary(m.out.all)$sum.all[,1:3]
## we’ll extract the matched dataset from the matchit object using match.data().
## m.data.all <- match.data(m.out.all)

list.matched <- m.out <- list()
for(j in 1:4){
    m.out[[j]] <- matchit(treatment.binary ~ male + college + age.young + conservative,
                          data = df.A1 %>% filter(gb==j),
                          replace = TRUE)  
    list.bal[[j + 1]] <- summary(m.out[[j]])$sum.all[,1:3]
    list.matched[[j]] <- match.data(m.out[[j]])
}
sink("balance_result.txt")
print(xtable(Reduce(rbind, list.bal))
      )
sink()


g1 <- love.plot(m.out[[1]], binary = "std", title="Pre-game")
g2 <- love.plot(m.out[[2]], binary = "std", title="Round 1 (draw)")
g3 <- love.plot(m.out[[3]], binary = "std", title="Round 2 (loss)")
g4 <- love.plot(m.out[[4]], binary = "std", title="Round 3 (win)")

pdf("result_balance_plot.pdf", family="sans",
    width     = 12, height    =10)
NetworkChange:::multiplot(g1, g2, g3, g4, cols=2)
dev.off()


################################################
## Redo the analysis using the matched data
################################################
df.matched <- Reduce(rbind, list.matched)
df.sub <- df.matched
ict.results.R0 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==1)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R1 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==2)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R2 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==3)), 
                           treat = "treatment.binary", J=5, method = "lm")
ict.results.R3 <- ictreg(y.binary ~ 1, data = data.frame(df.sub %>% filter(gb==4)), 
                         treat = "treatment.binary", J=5, method = "lm")
ict_model.R0 <- model.maker(ict.results.R0)
ict_model.R1 <- model.maker(ict.results.R1)
ict_model.R2 <- model.maker(ict.results.R2)
ict_model.R3 <- model.maker(ict.results.R3)


ict.all <- bind_rows(ict_model.R0 %>% mutate(model="Pre-game"),
                     ict_model.R1 %>% mutate(model="Round 1 (draw)"),
                     ict_model.R2 %>% mutate(model="Round 2 (loss)"),
                     ict_model.R3 %>% mutate(model="Round 3 (win)")) 


df.five <- ict.all %>% filter(term == "Outgroup") %>% mutate(term=model) %>%
    dplyr::select(term, estimate, std.error)
df.five$term <- c("Pre-game",  "Draw", "Loss", "Win")
df.five$term <- factor(df.five$term, levels = c("Pre-game", "Loss", "Draw", "Win"))



df.actual <- df.five 
df.actual$estimate.all <- c(NA, rep(df.five$estimate[1], 3))
df.actual$std.error.all <- c(NA, rep(df.five$std.error[1], 3))

g.matched <-  df.actual %>% 
    ggplot(aes(x=term, y= estimate, group=1)) + 
    xlab("")+ ylab("") +
    geom_point(size = 8, color="brown", alpha=0.5) +
    geom_pointrange(aes(x=term, ymin=estimate-1.66*std.error, ymax=estimate+1.66*std.error),
                    fatten = 0, color="brown", size=0.5, alpha=0.5) + 
    geom_line(size = 0.5, color="brown", linetype="dotted", alpha=0.7) +
    geom_point(data = df.five.null %>% filter(term == "Pre-game"), size = 8, color="brown", alpha=0.9) +
    ## average treatment
    geom_hline(yintercept= df.actual$estimate.all[2], alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] + 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    geom_hline(yintercept= df.actual$estimate.all[2] - 1.66*df.actual$std.error.all[2],
               linetype='dotted', alpha=0.3, color='brown', size=1) +
    annotate("text", x = 3, y = df.actual$estimate.all[2],
             label = "Average FIFA World Cup Effect") + 
    labs(title = "(b) Matched Data Result", y = "Outgroup Bias") + 
    theme(plot.title = element_text(face = "bold"),
          legend.position = "none")

pdf("result_matched_plot.pdf", family="sans",
    width     = 8, height    = 6)
g.matched
dev.off()
